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We study the contribution of advection by thermal velocity fluctuations to the effective 
diffusion coefficient in a mixture of two indistinguishable fluids. The steady-state diffusive 
flux in a finite system subject to a concentration gradient is enhanced because of long-range 
correlations between concentration fluctuations and fluctuations of the velocity parallel to the 
concentration gradient. The enhancement of the diffusive transport depends on the system 
size L and grows as ln(L/Lo) in quasi-two dimensional systems, while in three dimensions 
it grows as L^ 1 — L^ 1 , where Lq is a reference length. The predictions of a simple fluctu- 
ating hydrodynamics theory, closely related to second-order mode- mode coupling analysis, 
are compared to results from particle simulations and a finite-volume solver and excellent 
agreement is observed. We elucidate the direct connection to the long-time tail of the ve- 
locity autocorrelation function in finite systems, as well as finite-size corrections employed 
in molecular dynamics calculations. Our results conclusively demonstrate that the nonlinear 
advective terms need to be retained in the equations of fluctuating hydrodynamics when 
modeling transport in small-scale finite systems. 



I. INTRODUCTION 



Thermal fluctuations in non-equilibrium systems in which a constant (temperature, concentra- 
tion, velocity) gradient is imposed externally exhibit remarkable behavior compared to equilibrium 
systems. Most notably, external gradients can lead to enhancement of thermal fluctuations and 
to long-range correlations between fluctuations [TH3]. This phenomenon can be illustrated by con- 
sidering concentration fluctuations in an isothermal mixture of two miscible fluids, subjected to 
a macroscopic concentration gradient Vc. The solution of the linearized equations of fluctuating 
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hydrodynamics shows that concentration and density fluctuations exhibit long-range correlations, 
leading to a power-law divergence of the static structure factors for small wavenumbers k. When 
the species have different molecular masses and gravity g is present, the analysis predicts that 
fluctuations at wavenumbers below k g ~ (g ■ Vc) 1 ^ 4 are suppressed [6^11]. where c is the mass con- 
centration of the heavier species. Similar conclusions hold for fluctuations in a single-component 
fluid subject to a stabilizing temperature gradient [T2] . 

It is important to emphasize that this enhancement of large-scale (small wavenumber) concen- 
tration fluctuations occurs because of the non-equilibrium setting, and not because of the concen- 
tration gradient itself. Specifically, the enhancement is related to the dissipative flux through the 
system [13], which is zero at thermodynamic equilibrium. The top left panel of Fig. [T] shows a 
snapshot of the concentration field for a two-dimensional system in thermodynamic equilibrium in 
which there is a concentration gradient because of the sedimentation of the heavier species due 
to gravity, but no enhancement of the fluctuations. By comparison, the top right panel of the 
figure shows a non- equilibrium system with similar parameters but with an externally-imposed 
concentration gradient (via the top and bottom wall boundary conditions) and no gravity, reveal- 
ing much-enhanced fluctuations (noise) and large-scale features (clumping). If gravity is included 
in addition to the external gradient, the total diffusive flux is reduced and large-scale fluctuations 
(wavenumber k < k g ) are suppressed, as shown in the bottom left panel of Fig. [Tj In fact, if 
the same gravity as in the equilibrium case is imposed in addition to the external gradient, the 
total diffusive flux is essentially zero and the system is close to equilibrium again, giving no visible 
enhancement of the concentration fluctuations over the equilibrium shown in the bottom 

right panel of Fig. [T] These illustrative numerical results were obtained using a finite- volume solver 
for compressible fluctuating hydrodynamics |14j . 

The enhancement of concentration fluctuations is even more dramatic if the concentration gra- 
dient is at an interface, as in the study of the early stages of diffusive mixing between initially 
separated fluid components. As illustrated in Fig. [2j the interface between the fluids, instead of 
remaining flat, develops large-scale roughness that reaches a pronounced maximum until gravity 
or boundary effects intervene. These giant fluctuations [T2l [To] during free diffusive mixing 
have been observed using light scattering and shadowgraphy techniques [61 [SHU], finding good 
but imperfect agreement between the predictions of a simplified fluctuating hydrodynamic theory 
and experiments. In the absence of gravity, the density mismatch between the two fluids does not 
change the qualitative nature of the non-equilibrium fluctuations, and in this work we focus on the 
case of two dynamically-indistinguishable fluids. 
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Figure 1: Snapshot of the local concentration for a mixture of two ideal gases in a quasi- two-dimensional 
system, as obtained from a long run of a compressible finite- volume solver which includes thermal fluctuations 
[14j , with constant-temperature walls placed at the top and bottom boundary. The particles of species 1 (red 
end of the color scale) are four times heavier than particles of species 2 (blue end of the color scale) , and all 
systems are in mechanical equilibrium (i.e., the gravity force is balanced by the pressure gradient). (Top left 
panel) A system at equilibrium, in the presence of gravity go, with the heavier particles scdimcnting toward 
the bottom according to the equilibrium Gibbs-Boltzmann distribution. ( Top right panel) No gravity but 
a similar non- equilibrium concentration profile as in the top left panel is imposed via Dirichlet boundary 
conditions at the top and bottom walls, showing giant concentration fluctuations. (Bottom left panel) The 
same boundary conditions for the concentration are imposed as in the top right panel, but with gravity 
g = O.lc/oi showing a suppression of the large-scale giant fluctuations. (Bottom right panel) The same 
boundary conditions for the concentration are imposed as in the top right panel, but with the same gravity 
g = go as in the top left panel, showing no enhancement of the fluctuations. 

The giant fluctuation phenomenon arises because of the appearance of long-range correlations 
between concentration and velocity fluctuations in the presence of a concentration gradient. Based 
on nonlinear fluctuating hydrodynamic theory, it has been predicted that these correlations give rise 
to fluctuation-renormalized transport coefficients at larger scales |16H18j . However, the predicted 
contribution from fluctuations to transport at mesoscopic and macroscopic scales has only recently 
been computationally observed and reported by the authors [19|. This paper presents a detailed 
exposition of both the theoretical prediction for the enhancement of diffusion and the numerical 
simulations verifying these predictions. In particular, it is important to understand how the effective 



Figure 2: The rough diffusive interface between two miscible fluids at three points in time (top to bottom), 
starting from an initially perfectly flat interface (phase separated system), without gravity. Compare to the 
top right panel in Fig. [I] 

transport coefficients depend on the length scale of observation. This length scale may be related 
to the length Ax of the finite- volume cells used in a fluctuating hydrodynamic solver, or it may be 
related to the physical dimensions of a finite system such as a nano-channel transporting liquid or 
a nano-wire transporting heat. 

We consider diffusion in a mixture of dynamically identical but labeled (as components 1 and 
2) fluids [20] enclosed in a box of lengths L x x L y x L z , in the absence of gravity. Periodic 
boundary conditions are applied in the x (horizontal) and z (depth) directions, and impermeable 
constant-temperature walls are placed at the top and bottom boundaries. A concentration gradient 
Vc = {fiT—CB)/L y is imposed along the y axes by enforcing a constant mass concentration cy at the 
top wall and cb at the bottom wall. Because the fluids are indistinguishable, concentration does not 
affect the fluid properties, and the dynamics of the density, temperature and velocity fluctuations 
remains as in thermodynamic equilibrium. In this sense, concentration is passively transported by 
thermal fluctuations, analogous to diffusion of a passive tracer in a turbulent velocity field [21\ I22j. 
Note that for large L y the mass flux will be proportional to the self-diffusion coefficient of a tagged 
particle, independent of the magnitude of the gradient [20] . 

Since species are not changed in particle collisions, the diffusive transport of particle label 
(concentration) can only occur via advective motion of the particles between collisions. Kinetic 
theory shows that at steady state the particles of a given species (denoted either with a subscript 
or with a parenthesis superscript) have a non-zero macroscopic momentum density j\ = Pl^i 
j = —P2V2, where p denotes density and v velocity. If the labeling of the species is ignored, the 
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system is at equilibrium and the overall center-of-mass velocity vanishes, v = cv\ + (1 — c) t>2 = 0. 
More detailed kinetic theory [23|, 123] shows that the inter-species velocity v\2 = v\ — v<i quickly 
relaxes to its equilibrium value, 

giving the Fickian diffusive flux for the mass concentration c = p\j p [23J, 

3 1 = P\V\ = p\v - px (Vc) = -px (Vc) , 

where x is the mass diffusion coefficient. The local fluctuations around the macroscopic mean, 
Pi = Pi + 5 pi and v\ = v\ + 5vi, can also make a non-trivial contribution to the average mass flux, 

til) = (Pivi) = -PX (Vc) + ((S Pl ) (5 Vl )) , (2) 

if they are correlated, which in fact they are in the presence of a concentration gradient. 

The fluctuating hydrodynamics formalism [5J, [25] is the most direct way to calculate steady-state 
correlations between hydrodynamic variables, especially if a spatial Fourier transform is used to 
separate different wavenumbers. Converting the correlations from Fourier to real space requires 
integrating over all wavenumbers, which gives qualitatively different results in two and in three 
dimensions, as detailed in Section [ITJ In two dimensions the average mass flux, or effective diffu- 
sion coefficient, is found to grow logarithmically with system size, while in three dimensions an 
asymptotic "macroscopic" value is reached for sufficiently large systems. 

As seen from Q, the correlation that is needed is that between the density and velocity of a 
given species, dpi and 5v\. However, in the standard "single-fluid" hydrodynamic description of 
mixtures, unlike the little-understood "two-fluid" models [23l [23] , the individual species velocity V\ 
(or equivalently, V\2j is not maintained as an independent variable and instead only the center-of- 



mass velocity v appears [5]. In Section III we report results from simulations based on the Direct 



Simulation Monte Carlo (DSMC) particle method [26, 27J, in which we calculate the spectral 



correlations between 5p\ and 5v\ and also the mass flux. The results presented in Section III A 
indicate that the correlation between 5p\ and Svi is well- approximated by the prediction of the 
incompressible single- fluid theory presented in Section [TTl The effective diffusion coefficient is found 



to increase with system size in accordance with the theory as well, as detailed in Section III B For 
systems with aspect ratio close to unity the use of the periodic (Fourier-based) theory is not 
appropriate and the proper boundary conditions ought to be taken into account, as we do by using 



a recently- developed compressible finite- volume solver |14j in Section III B 1 Very good agreement 
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is observed between the finite-volume simulations and the particle results over a broad range 
of system sizes, once the local diffusion coefficient that appears in the equations of fluctuating 
hydrodynamics is adjusted to match particle data for a chosen reference system. This locally- 
renormalized diffusion coefficient is also measured in the particle simulations and found to match 
the theoretical predictions reasonably well, as discussed in Sec. |III B 2 



The present paper builds on an extensive prior literature on the renormalization of the diffusion 



coefficient by hydrodynamic fluctuations and interactions. In Section IV we discuss connections 
to prior work in more detail, and find that several different theoretical approaches produce the 
same results as the very simple, intuitive, yet extensible fluctuating hydrodynamic theory. In 
particular, we compare to previous mode-mode coupling theories of the long-time tails in the 
velocity autocorrelation function, as well as to theories of finite-size effects on diffusion in periodic 



systems. Furthermore, in Section IV A we re-examine existing data from hard-disk molecular 
dynamics simulations to find that the simple theory describes the system size dependence of the 
long-time diffusion coefficient of hard disks as well, confirming that the phenomenon we study is a 
generic property of particle fluids and not an artifact of DSMC. For large system sizes, however, 
we find that a more sophisticated self-consistent theory is necessary, and make some preliminary 
attempts at an explicit self-consistent calculation in both two and three dimensions, before offering 
concluding remarks and discussing future research directions. 



II. FLUCTUATING HYDRODYNAMICS 



At mesoscopic scales the hydrodynamic behavior of fluids can be described with continuum 
stochastic PDEs of the Langevin type [28, 29J. Thermal fluctuations enter as random forcing terms 
in the Landau-Lifshitz Navier-Stokes (LLNS) equations of fluctuating hydrodynamics [25\ !3Uj . 
For a mixture of two indistinguishable fluids, neglecting viscous heating, the compressible LLNS 
equations are I3"T] 

D t p = - p(y -v) 

p (D t v) = — VP + V • (r]Vv + E) 
pc v (D t T) = - P (V • v) + V • (kVT + H) 

p(Ac)=V-[px(Vc) + *], (3) 

where D t O = 9 t D + v ■ V (□) is the advective derivative, V« = (Vi> + Vu T ) — 2 (V -v) 1/3, and 
the pressure is P = p(kBT/m) = pc^, where ct is the isothermal speed of sound. The viscosity 
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rj, thermal conductivity k, and the mass diffusion coefficient x may in general depend on the 
state. The capital Greek letters denote stochastic fluxes that are modeled as white-noise random 
Gaussian tensor and vector fields, with amplitudes determined from the fluctuation-dissipation 
balance principle [32J, 

S = v^tpW, * = y/2mxp c(l - c) W. and S = T y / 2 K k B VV (4) 

where m is the fluid particle mass, and W, W and VV are white-noise random Gaussian tensor 
and vector fields with covariance 

(Wij(r, t)V\&(r', 0) = (fofy + *a*i* - 2<%<W3) - 0*(»" " r ')> 
(vV i (r,t)VV*(r',t')) = (>^i(r,t)V^(r',0) = (<%) *(i - i')<5(r - r'), 

where star denotes complex conjugate. 

In addition to the usual Fickian contribution, the flux in the equation ^ for c = c + 5c includes 
advection by the fluctuating velocities, v = v + 5v = 5v. Ignoring density fluctuations, 

d t (5c) + (5v) ■ (Vc) = V • [ X V (5c) - (5c) (5v)] + p' 1 (V • *) . (5) 

Interpreting the non-linear stochastic partial differential equation (SPDE) ^ requires some form 
of regularization (smoothing) of the stochastic forcing, usually approached using a perturbative 
approach |16H18|, [33] . To leading order, we can approximate the advective contribution to the aver- 
age diffusive mass flux, using the solution of the linearized equations of fluctuating hydrodynamics, 
which can be given a precise meaning |34j . Specifically, we anticipate a relation of the form 

- ((5c) (5v)} « - ((5c) (5v)) hncav = (A X ) Vc, 

leading to an effective diffusion coefficient Xeff = X + that includes an enhancement Ax due 
to thermal velocity fluctuations, in addition to the bare diffusion coefficient x- I n Appendix [A| we 
give some simple estimates of the relative magnitude of in relation to x, demonstrating that 
the enhancement due to velocity fluctuations is expected to be much larger for dense liquids than 
for dilute gases. 

A. Fluctuation-Enhanced Diffusion Coefficient 

In order to analyze the stationary solution to the linearized equations of fluctuating hydro- 
dynamics, we will apply a Fourier transform in all directions as done in Ref. [18], even though 
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the direction of the gradient is not periodic. One can justify this approximation by considering a 
periodic background concentration field, maintained at steady state via some external potential, 
and then calculate the mass flux in the vicinity of the plane y = L y /2 as a function of the local 
concentration gradient in the limit of infinite period L y [5]. 

To simplify the analysis, we can neglect density and temperature variations, p = po and T = Tq, 
to obtain the isothermal incompressible approximation, 



d t v =V [-v ■ Vv + uV 2 v + p' 1 (V • S) 
d t c = - v • Vc + xV 2 c + p' 1 (V • *) . 



(6) 
(7) 



dB, 



(8) 



where v = rj/p, and v • Vc = V • (cv) and v • Vt? = V • (tm r ) because of incompressibility, 
V • v = 0. Here "P is the orthogonal projection onto the space of divergence-free velocity fields, 
"P = I — k~ 2 (kk*) in Fourier space (denoted with a hat). 

The linearized form of (6j7) in the Fourier domain is a collection of stochastic differential equa- 
tions, one system of linear additive-noise equations per wavenumber, of the form 

1/2 

2p- 1 uk B Tk 2 V 

2p~ 1 x?7ic(l — c) k 2 

where 3 is a collection of independent Wiener processes. At steady state the correlations between 
the Gaussian fluctuations are described by the matrix of static structure factors (covariance matrix) 

'(5v)(5v)*\ ((fo)(6c 
'(Sc)(Siy) ((6c)(Sc 

The static structure factor matrix consists of a short-ranged equilibrium contribution and a long 
range non-equilibrium contribution, 





5c 




uk 2 V 




5c 




d 


6v 




9c Xk 2 




5v 


dt + 



+ Vc 







AS C , V (AS c , c )Vc 



p- l k B TV 

mp~ l c(l — c) 

The explicit form of S can be obtained as the solution of a linear system derived from (|8j) using 
the stationarity condition dS = [Hj. The concentration fluctuations are enhanced as the square 
of the applied gradient [18] , 



(5c)(5cY 



2 _ ( S [ n 2 0: IT-...\2 



(A5 C)C ) (Vc) 



e) (vc 



(9) 



while the correlation between the concentration fluctuations and the fluctuations of velocity parallel 
to the concentration gradient are linear in the applied gradient |18j . 

k B T 



(5c)(5v 



AS r . 



Vc 



p{u + x)k 2 



(sin 2 (9) Vc. 



(10) 
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where 9 is the angle between k and Vc, sin 2 9 = k\jk 2 . The power-law divergence for small 
k indicates long ranged correlations between 5c and 6v and is the cause of both the giant fluc- 
tuation phenomenon and the diffusion enhancement. As seen from Q, the actual correlation 
that determines the diffusion enhancement is S n ji) = /(5c)(5v^ )*\, which is approximated as 



S n „(i) ~ pSc,v« in (p 



Ill A 



7); this approximation is discussed and justified in Section 
The mass flux due to advection by the fluctuating velocities can be approximated as 

(5c (r, t) 6v (r, t)> = (2vr)^ 6 J dk dk' (si: (k, t) 5v* (k' , t) } e < k - k '> r = (2vr)" 3 J 5 C ,„ (fc) dk, 

(11) 



which together with ( 10 ) gives an estimate of the diffusion enhancement 

A X = - (2vr)- 3 f AS C , V ,, (k) dk = [ ( sin 2 9) AT 2 dk. (12) 

Because of the £; _2 -like behavior, the integral over all k above diverges unless one imposes a lower 
bound, fcmin ~ 27r/L in the absence of gravity, and a phenomenological cutoff k max ~ ir/L mo i [18] 
for the upper bound, where L mo \ is an ad-hoc "molecular" length scale. Importantly, the fluctuation 
enhancement Ax depends on the system size L because of the small wavenumber cutoff. 

1. Two Dimensions 

For a quasi two-dimensional system, L z <C L x <C L y , we can replace the integral over k z with 
2tt/L z and integrate over all k y . This leads to an average total mass flux that grows logarithmically 
with the system width L x for a fixed height L y , 

ksT r\k x \<TT/L mo i £ 2 ksT L x 

A X (L X )^ {2 ^2 p{x + u)Lz J lkxl > 2n/Lx ( k 2 + k 2f dkydkx = 4irp(x + v)L z ^ 2L^{ (13) 

When the system width (perpendicular to the gradient) becomes comparable to the height (parallel 
to the gradient), boundaries will intervene and for L x 3> L y the effective diffusion coefficient 
must become a constant, which is predicted to be a logarithmically-growing function of L y in two 
dimensions. 

It is important to emphasize here that the chosen value of L mo \ is arbitrary. The hydrodynamic 
theory models the effective diffusion coefficient as the sum of the "bare" diffusion coefficient x 
and the "enhancement" A%, but the two cannot be separated because every measurement must 
be performed for some finite L x . One can thus simply define x to be the value of the measured 
diffusion coefficient for some reference width Lq > 2L mo \, and predict that for L x > Lq, 

J 2D ) ~ y + k JL i n h ( U ) 
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For this prediction to be accurate, however, Lq ought to be chosen to not be too large, so that 
the enhancement of the diffusion relative to the "molecular" contributions is small and simple 
quasi- linearized theory applies, but also not too small so that fluctuating hydrodynamics applies. 

Because we are explicitly concerned with the effect of a finite width L x , the integral over k x 
should be replaced by a discrete sum over the wavenumbers consistent with periodicity, k x = 
k x ■ 2ir/L x , where k x 6 Z. If one calculates the difference between a system of width 2L X and a 



system of width L x , then it is easily seen that the integral over k x in Eq. (13) ought to be replaced 
with the following sum over n x , 



(27^ 



f(k x ) dk x - / f(k x ) dk x 

\k x \>n/L x J\k x \>2ir/L x 



<— > L- 1 £ [f{2n x -l)-f{2K x )\ 

K x y£0 



Even though f (k x ) ~ k^ 1 in (13) is not integrable, the difference in the square bracket above 
goes like k~ 2 and the sum can be done explicitly, giving exactly the same answer as the integral 
estimate, 

2. Three Dimensions 

Now we consider a system where L x = L z = L L y , and study how the effective diffusion 
coefficient changes with L. In three dimensions, the relative contribution from large wavenumbers, 
i.e., small scales, is larger than in two dimensions. We can use the integral approximation to 
examine the asymptotic behavior for large L, 

^ k B T r\^\<^ k 2 + k 2 ^ ^ ^ = In (1 + V2) k B T / 1 _ 1 

( 271 ") 3 P(X + ^) J\k x/z \>2Tv/L {kl + kl + klf V X " MX + V) V2^mol L 

We see that in three dimensions Xeff converges as L — > oo to the macroscopic diffusion coefficient, 
but for a finite system the effective diffusion coefficient is reduced by an amount ~ L _1 due to the 
truncation of the velocity fluctuations by the confining walls, 

(3D) ak B T (15) 

Calculating the exact value of a requires performing a sum over k x and k z instead of integrals over 
k x and k z , as we have done numerically. The numerical results suggest that, as in two dimensions, 



the difference in Xcs*^ between two systems attains a finite value as L mo \ — > 0, justifying (15) for 
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III. PARTICLE SIMULATIONS 

This sections verifies the predictions of fluctuating hydrodynamics by particle simulations. Here 
we employ the Direct Simulation Monte Carlo (DSMC) particle algorithm |26t I27|. in which de- 
terministic interactions between the particles are replaced with stochastic collisions exchanging 
momentum and energy between nearby particles. The collision rules ensure local energy and mo- 
mentum conservation and a thermodynamically-consistent fluctuation spectrum [65\ I36|. Previous 
careful measurements of transport coefficients in DSMC using nonequilibrium methods have been 
limited to quasi one-dimensional simulations, in which there is only one collision cell along the 
dimensions perpendicular to the gradient [37 . The effect we are exploring here does not appear in 
one dimension as it arises because of the presence of vortical modes in the fluctuating velocities. 

We have performed DSMC calculations for an ideal hard-sphere gas with molecular diameter 
a = 1 and molecular mass m = 1, at an equilibrium density of po = 0.06, with the temperature 
kept at ksTo = Tq = 1 via thermal collisions with the top and bottom walls. A uniform concen- 
tration gradient along the vertical (y) direction is implemented by randomly selecting the species 
of particles to be one with probability ct/b when they collide with the top/bottom wall. Each 
DSMC particle represents a single hard sphere so the mean free path is A = 3.75 and the mean 
free collision time is r = 2.35. The DSMC time step was chosen to be At = r/2, and the collision 
cell size is either Ax c = A or Ax c = 2A. 

The DSMC algorithm simulates a dilute gas, for which the enhancement of diffusion is weaker 
than for dense fluids (see Appendix [A]) . Nevertheless, the computational efficiency of DSMC makes 
it preferable to molecular dynamics for this study. The DSMC method employed here uses a grid 
of collision cells, thus introducing discretization artifacts into the particle dynamics. While it is 
possible to eliminate these grid effects entirely \36\ [38] , the associated increase in computational 
cost and the difficulty of parallelization would make some of the large-scale particle simulations 
presented here infeasible. Furthermore, decreasing the density in order to increase the mean free 
path and reduce the grid effects would make the relative size of the effect we are trying to observe 
too small compared to statistical errors. We have verified that quantitatively identical results 
are obtained for two different choices of DSMC collision cells, Ax c = A or Ax c = 2A, once the 
discretization correction to Chapman-Enskog kinetic theory for the transport coefficients is taken 
into account [39 41j . 

In addition to the DSMC collision cells, which determine the microscopic dynamics of our 
particle simulations, obtaining hydrodynamic quantities such as velocity requires using a grid of 
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N x x N y x N z sampling or hydrodynamic cells, each of volume AV = |AV| = Ax Ay Az. The 
sampling of hydrodynamic quantities is performed every 5 DSMC time steps, at a snapshot time 
that is randomly chosen. Sampling at random time intervals ensures that there is no measurement 
bias due to the lack of time invariance in the particle dynamics, and gives similar results as sampling 
at the mid point of each time step [37]. At each snapshot, we obtain the instantaneous mass ttjav = 
(Ay) pav an d momentum Pav = (AV) j av in each sampling cells by adding the contributions 
from all particles inside the given sampling cell. We can do this sampling taking into account 
either all of the particles, pav an d Jav> or j us ^ particles of the first species, which we indicate by 
a species subscript or parenthesis superscript, p^y and j'av- ^ or eacn sampling cell, we obtain an 
instantaneous velocity i>av = Jav/pav (similarly for v^L) and mass concentration c = p^v/PAV- 
We obtain discrete static structure factors (spectral correlations) from time averages of products 
of discrete Fourier transforms (DFTs) of the instantaneous variables. For comparison between the 
particle simulations and the theory we use a reference length Lq = 16A. 



A. Static Structure Factor 



In order to compare the prediction (10) to results from particle simulations, we need to convert 



the continuum static structure factor S CjV ,, (k) into a discrete structure factor S c>v ., (re) for finite- 
volume averages of the continuum fields. Here the set of N x x N y x N z wavenumbers k £ Z 3 indexes 
the discrete set of wavevectors compatible with periodicity, k (k) = 2ir (re^L" 1 , n y L~ x , k z L~ x ^. A 
relatively straightforward calculation shows that 

S ^w (^) = F ^ i k i k K)] . ( 16 ) 



where the sum goes over all resonance modes, k! = (k x + N x Ak x , K y + N y An y , k z + N z Ak z ) for 



all Ak g Z 3 , and F AV (k) = F x (k x ) F y (k y ) F z (k z ) is a product of low-pass filters of the form 



F x (k x ) = sine 2 (k x Ax/2) , (17) 



where sinc(x) = sin(x)/x. The sum in (16) can easily be evaluated numerically because the terms 



decay rapidly [c.f. (10)] 



In Fig. [3] we compare the theoretical prediction for pS CjV ,, (n) to results from particle simulations 
for the discrete structure factor 



S m(«) - ( (<Vj ) ( Sv 1} 
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A k y /(2ir) 




\k y /(2ir) 



Figure 3: (Top) Discrete structure factor S (i) from quasi two-dimensional DSMC runs with L x = 64A, 

P1 ' v i\ 

L y = 512A and L z — 2A, for several wavenumbers k x — k x ■ 2tt/L x (circles k x — 2, squares n x = 6, diamonds 
k x = 10, triangles n x = 14), compared to pS c ^ v ^ (solid lines of the same color) as predicted by the linearized 
periodic theory ( 10|16 1. The sides of the DSMC collision cells are Ax c = Ay c = 2A = 7.5. Note that for a 
fixed k x we expect the structure factor to decay as k~ 4 . (Bottom) All the parameters, including the system 
and sampling cell size, are as for the system in the top panel, but now the DSMC cells are twice smaller, 
Ax c = Ay c = A = 3.75. Note that this change of the DSMC cell size changes the kinetic theory prediction 
for viscosity by more than 20% [55] . 
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for two different sizes of the DSMC collision cells. The fact that there is little difference between 
the two panels in the figure verifies that the details of the microscopic collision dynamics do not 
affect the mesoscopic hydrodynamic behavior. In Fig. [4] we plot the discrete structure factor from 
the particle simulations for wavevectors perpendicular to the gradient (i.e., k y = 0), for systems of 
different width L x . 




Figure 4: Discrete structure factor v (i) from quasi two-dimensional DSMC runs with L y = 512A and 
L z — Ax c — Ay c — 2 A, for wavevectors perpendicular to the gradient {k y = 0). Results for systems of 
different width L x are shown with symbols (see legend). For comparison, the theoretical prediction for 



pS c>v ^ for an infinite periodic system (10 16 ) is shown with a line. Deviations from the predicted k x 2 
power-law divergence are clear for k x < 2ir/L y w 3 • 10 -3 due to the influence of the top and bottom walls. 



It is expected that compressibility effects would affect ^(i). Indeed, this is what we observe 
in simulations, however, as Fig. |3| demonstrates, the incompressible isothermal theory for pS CjV ^ 
is in very good agreement with particle data for . Good agreement between the simulation 

data and the simple theory is also seen in Fig. [4j except for k x comparable to 2n/L y . As expected, 
for the smallest wavenumbers the top and bottom walls intervene and the actual correlation is 
smaller than the predicted k~ 2 divergence. 

In order to construct a theoretical prediction for ^(i), one must not only include the effects 
of compressibility but also replace the "one-fluid" approximation (|3]) with a corresponding "two- 
fluid" compressible hydrodynamic theory |23} [24] . This can be seen by noting that the fluctuating 
equations ([3]) assume that relation Q applies to the fluctuating V\2 an d c instead of their means. 
Such an assumption leads to unphysical bias of order ((5c) 2 ) in the mean inter-species velocity 



15 



i v i2) because of the nonlinearity in the denominator c(l — c). In fact, the fluctuations 5v\2 
and 5c should be uncorrelated, as seen from a two-fluid fluctuating theory. Here we use the 
incompressible isothermal approximation for pS c , v ^ as a proxy for S ^ in order to construct 
theoretical predictions for the diffusion enhancement. 

B. Fluctuation-Enhanced Diffusion Coefficient 

As we already explained, instantaneous hydrodynamic quantities, denoted with a subscript 
AV, are sampled from the particle data by taking snapshots of the particle state using a grid of 
sampling cells of volume A V . The ensemble average of a given quantity, which we will denote with 
angle brackets, is obtained by averaging over many snapshots once a steady state is reached, and 
additional averaging can be performed over all sampling cells with the same y position since the 
steady state averages cannot depend on x and z. We estimate the macroscopic mean mass density 
p, partial density p%, partial momentum density j 1 , partial velocity V\ and concentration c as 

P = Po = (pav) , Pi = (pav) = d P and ix = jo } = (j'av) = P 1 ^ 1 - 

We also define the mesoscopic velocity and concentrations to be the ensemble averages of the 
instantaneous values, 

v ^ = ( v av) and c o = ( c av) , 

where the subscript zero will be used to simplify the cumbersome notation. It is important to point 
out that for non-conserved quantities such as v and c the mesoscopic mean can be different from 
the macroscopic mean due to fluctuations [32l S3], v\ ^ and c 7^ co. For conserved quantities 
(e.g., ji and jfg 1 ), however, the mesoscopic and macroscopic ensemble means are equal and in fact 
independent of Ax and Az (but not necessarily Ay). 

In particle simulations, we calculate the effective diffusion coefficient XeS from the momentum 
density of one of the species along the vertical direction, 

J| = POXett 7 T— ~ POXett Vc, (18) 

II Ly - Ay 

where we measure ct and cb in the top and bottom layer of sampling cells (whose centers are a 
distance L y — Ay from each other) to empirically account for the small concentration slip in DSMC 
(about 0.5% with these parameters). Numerical experiments have verified that j^p matches the flux 
obtained from counting the average number of color flips at the top or bottom walls. Furthermore, 
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the results verify that Xeff is essentially independent of the magnitude of the concentration gradient, 
and that the change in the effective gradient Vc as L x or L z is changed, keeping L y fixed, is much 
smaller than the change in Xeff- 

The traditional definition of a "renormalized" diffusion coefficient [16\ [T7] as the macroscopic 
limit of XeS, only works in three dimensions and is not very useful for confined systems. Instead, 
for each sampling cell, we define a locally renormalized diffusion coefficient xo y i & 

,W„(1) _ / JXi \ /,.(!) \ 



P*"°o' = \Pav) { v av) = PXO (Vc) , (19) 

where we have accounted for the fact that the macroscopic concentration gradient dc/dy may 
depend on y. In fact, such a dependence is observed in the particle simulations, and we have 
approximated the local concentration gradient dc/dy by a numerical derivative of a polynomial 
fit of degree five to c(y). Figure [5] shows that the empirical xo is independent of y, except for a 



boundary layer close to the top and bottom walls. This is an important finding, since (19) is a 



constitutive model that is assumed to hold not just at the macroscale but also at the mesoscale, 
notably, xo is an input parameter for fluctuating hydrodynamics finite- volume solvers |14j . 
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Figure 5: The renormalized diffusion coefficient xo, as defined by Eq. (19), as a function of the y position 



of the sampling cell for DSMC systems with several system widths L x . The estimate of xo shown in Fig. [6] 
is simply the average over all sampling cells further than 10Ay away from the top and bottom walls. 



Figure [6^, shows how the effective Xeff an d renormalized xo diffusion coefficients change as the 
width of the system L x is increased while keeping the height L y fixed for two different quasi two- 
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dimensional DSMC systems. For System A, the DSMC collision cells are cubes of side Ax c = 7.5 = 
2A, while each sampling cell contains 2x2x1 collision cells, or N p = 101 particles on average. The 
height of the box is L y = 256A = 960 and the imposed concentrations at the walls are cb = 0.25 
and ct = 0.75. For System B, the DSMC parameters and ct/b are the same as for to System A, 
but the sampling cells are twice as large, 4x4x1 collision cells each, and the system height is 
twice as large, L y = 512A = 1920. We obtain similar results using twice smaller collision cells (not 
shown). For the quasi two-dimensional systems, the thickness is L z = 7.5 = 2A and there is only 
one DSMC collision cell along the z direction. Figure [6^, shows that Xefi grows like lnL x , with 



a slope that is well-predicted by Eq. (14). For widths larger than about 8 mean free paths, xo 
becomes constant and rather similar to the kinetic theory prediction. It is important to point out 
that xo is n °t a fundamental material constant and in fact depends on the shape of the sampling 



cells (see Section III B 2 ) . 

In Fig. [6}) we show results from three dimensional DSMC simulations, in which the system 
width (x) and depth (z) directions are equivalent, L z = L x = L, and the rest of the parameters 
are the same as for System A. Similar behavior is seen as in two dimensions, except that now 
the effective diffusion grows as — L _1 and saturates to a constant value for large L, assuming that 
L y > L. 



1. Corrections due to finite height 



The predictions of the simplified fluctuating hydrodynamic theory, Eqs. ( 14 ) and ( 15 ) , are shown 



in Fig. [6] and seen to be in very good agreement with the particle simulations for intermediate 
L x . However, the particle data shown in Fig. [6^, shows measurable deviations from the simple 
theory for L x > L y /2. To understand the discrepancy, recall that the incompressible isothermal 
theory assumed that L y is essentially infinite and thus in two dimensions XeS grows unbounded in 



the macroscopic limit. A scaling analysis suggests a modification of (14) to account for the finite 
height of the system, 



y (2D) 



x + 



k B T 



ln~ - f(r) 
Lo 



(20) 



4irp(x + v)L z 

where r = L x /L y is the aspect ratio of the system, and f(r) is some function that is close to zero for 
small r and grows asymptotically as ln(r). Therefore, when L x 3> L y , Xeff saturates to a constant 
value that grows as ln(L y /Lo)- 

One can extend the theoretical calculations to account for the hard wall boundary conditions 
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Figure 6: (Panel a, top) The effective \eS and the renormalized xo diffusion coefficients as a function of the 
width of the system L x in two dimensions. Numerical results for System A (DSMC and SPDE) and System 
B (DSMC) are shown with symbols (see legend). The error bars for all of the numerical data are comparable 



or smaller than the size of the symbols. The theoretical predictions (23 1 are evaluated numerically and 
shown with lines. The bare diffusion in the theory and SPDE calculations is adjusted so that for L x = 16A 
the effective diffusion is the same as that measured in the particle simulations. (Panel b, bottom) Same as 
the top panel but in three dimensions. The inset highlights the behavior. 
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in the y direction [5], however, such a calculation is non trivial. Instead, we have used the finite- 
volume solver developed in Ref. [H] to solve the non-linear system of SPDEs Q for the same 
system dimensions as in the particle simulations. The results, shown in Fig. |6| are in excellent 
agreement with the particle simulations for the larger system sizes. Note that while our SPDE 
solver includes all of the nonlinear terms in ([3]), we may artificially reduce the amplitude of the 
noise and thus the magnitude of the fluctuations by some factor eCl, This reduces the effect of 
the nonlinearities and effectively gives a quasi-linearized finite- volume SPDE solver. The advective 
mass flux due to the velocity fluctuations can be estimated as 

Ajj 1 ^) « 6- 2 ((6 P1 ) « e-*p((Sc) (&,„)) , (21) 

and may depend on y especially close to the walls or when L y > L x . The sum of the average 
diffusive and advective mass fluxes must be independent of y, 

^(i) _ c T -c B _ dc{y) -.(i) 

jj -PXeff— ^ — = PXo-^ - ^ (y), (22) 
which implies that the macroscopic concentration profile c(y) is affected by the fluctuations as well 



and cannot be strictly linear. From the conditions c(0) = cb and c(L y ) = ct and (22) we obtain 
the relation 



XcS 



Jv=0 L 11 



dy, 



>y=o 

which is how we calculate the effective diffusion coefficient from the numerical SPDE solution. We 
have verified that the results are independent of e to within statistical accuracy for e < 1/2. 

The velocity-concentration correlation Ai| (y) obtained from the finite- volume solver is shown 
in Fig. [7J along with the corresponding particle data for comparison. Excellent agreement is seen, 
demonstrating that the finite- volume solution correctly accounts for the influence of the boundaries. 
Note that the partial velocity V\ is not included as an independent variable in ([3j and the mean 
velocity v is used instead. When compressibility effects are included, v, unlike Vx, is correlated 
with pi even in one dimension. This makes a direct comparison between the effective diffusion 
coefficient in the particle and finite- volume simulations difficult. However, the dependence of Xeff 
on system size should be the same in both types of simulations, once the bare transport coefficient 
is adjusted empirically. 



2. Renormalized Diffusion Coefficient 



The renormalized diffusion coefficient \o m the Fickian diffusive flux is an input to the SPDE 
calculations and assumed constant. In our calculations we used the prediction of kinetic theory 
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Figure 7: The fluctuating contribution to the mean diffusive flux, ^(<5pi) f^f^My, as a function of the 
height of the sampling cell y, for several system widths L x , keeping L y — 256A. Data from quasi-two 
dimensional particle simulations (System A) is shown with symbols, and lines of the same color show data 
for p ((<5c) (<5«||)) obtained using a two-dimensional finite- volume solver P3| for the LLNS equations |3j). The 
error bars for the particle data are shown for L x = 4A, and are similar or smaller for the remaining systems 
and for the finite- volume results. 



also shown in Fig. |6| In finite- volume solvers, the spacing of the computational grid plays 
the equivalent of the cutoff length L mo \, and therefore the effective mass flux depends on the grid 
spacing. Furthermore, there are numerical grid artifacts in the SPDE solution at length and time 
scales comparable to the numerical discretization parameters |14j . To correct for these errors, we 
have added a constant to the effective diffusion coefficient obtained from SPDE runs to match Xeff 
from the particle simulations for L x = Lq = 16A. This correction essentially renormalizes xo based 
on the size of the finite- volume hydrodynamic cells. 



One can think of xo defined via ( 19 ) as the fluctuation-renormalized diffusion coefficient at length 



scales determined by the shape of the sampling or observation volume AV. In this sense, xo is the 
physical-space equivalent of the wavenumber-dependent diffusion coefficient x (k,cu = 0) commonly 
used in linear response theories |1 6|. I20j . A theoretical prediction for xo can be obtained by starting 
from linearized theory for the fluctuating fields p\ (r,t) = p\ + bp\ (r,t) and v\ (r,t) = v\ + 
5v\ (r, t). The instantaneous velocity in a given sampling cell was defined through the instantaneous 
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momentum density j 1 = p\V\ averaged over the sampling cell, 



o {1) 
Pav 



V AV /-Pi \ V AV / 



Pi 



Jav Pi v i dr 
Iav P 1 dr 



Jl) ,,(1) \ A,' 



where to second order in the fluctuations, 



Aj F = AV~ 2 dr / dr' (pi{r,t)vi(r',t) 
Jav Jav 



--AV- 



(2tt)- 

(27T)- 



'AV 
r/r 



AV 



dr' (2^r 6 / dfc 



AV 



AV- 



dr 



AV 



k' 

dr' e ik < r - r '^ 



AV 



dk' \5pi (k, t) Svi 
"Spi,vi (k) dk 



(k',t) 



i(k-r—k'-r') 



' I F AV (k) S pi>Vl (k) dk, 

Jk 



and Fav (k) is the low pass filter that already appeared in Eq. (17). The result of this calculation 



[c.f. Eq. (12)], 



Xeff = X 

Xo 



(27T) 



p^AS 



(1) 



(k) 



dk 



X - (27T)- 3 f[l- F AV (fc)] Ipo'AS (1) 
Jk L p1 ' ii 



(fc) 



dfc, 



(23) 



shows that Xeff includes contributions from all wavenumbers present in the system, while xo only 
includes "sub-grid" contributions, from wavenumbers larger than 2ir/ Ax. The theoretical predic- 
tions shown in Fig. [6] are based on numerically evaluating ( 23 ) after replacing the integrals over 
k x and k z (23) with the appropriate sums, assuming S pijVl ~ pS c , v and using Eq. (10). The bare 
diffusion coefficient x is adjusted so that Xeff matches the particle result for L x = Lq = 16A, and 



good agreement is observed between (23) and the particle data for xo f° r ah but the smallest L x . 

While it is intuitive to expect that the bare diffusion coefficient should account for molecular, 
or non-hydrodynamic, degrees of freedom, the division XcS = X + ^x is arbitrary, and in fact there 
is no unambiguous way to define x- This is evident in the theory because of the need to introduce 
an ad-hoc molecular cutoff as a way to separate the "microscopic" from the "mesoscopic" scales. By 



contrast, the locally renormalized diffusion coefficient xo defined in (19) explicitly depends on the 
size of the sampling (hydrodynamic) cells AV = \AV\ used to define the hydrodynamic quantities 



from the particle configuration. Combining the two equations in (23) gives the renormalization 
relation at large scales, 



Xeff = Xo (AV) - (2tt)- 



Fav (k) 



p^AS (1) (fc) 



dk, 



which eliminates the dependence on the ad-hoc cutoff wavenumber since Fav filters contributions 
from large wavenumbers, at least within the simple perturbative (quasi-linear) theory. 
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IV. CONNECTIONS TO EARLIER WORK 



While our computer simulations are the first hydrodynamic study of the dependence of trans- 
port on system size, there is a substantial body of literature that has discussed the effect from a 
theoretical perspective or studied smaller particle systems. In this section we explicitly connect 
our analysis to previous approaches, and discover direct relations with work that might have, at 
first sight, been assumed to be unrelated. 



A. Relation to Long-Time Tails 

It is well known that the self- diffusion coefficient is given by the integral of the equilibrium 
velocity autocorrelation function (VACF) C{t) of the fluid particles [3j. The long-time tail of 
C(t) has been extensively studied in the literature both computationally and through several 
theories, including heuristic hydrodynamic arguments [444 I45| . kinetic theory [IB] and (second- 
order) mode-mode coupling hydrodynamic theory |47j . Ultimately all derivations give the same 
result including not just the power-law dependence but also the coefficient of the tail, specifically, 
in three dimensions C(t) ~ k B T / |l2p [ir (x + v ) ^] 3 ^ 2 |i while for quasi-two dimensional systems 
C(t)^k B T/ [8tt P L z ix + v)t\. 

A crucial point is that the VACF explicitly depends on the system size due to periodic 
boundaries, and so its integral, which gives the diffusion coefficient, also depends on system 
size. More explicitly, ignoring acoustic effects, the VACF has the power-law dependence only 
for L^ ol / (x + v) <C t -C L 2 / (x + v )i an d it decays exponentially for large times [38]. Ignoring 
prefactors, the contribution of the tail to the diffusion coefficient in three dimensions is estimated 
as 

■"4ol/(x+«0 P X + V) t 3/2 P (X + ") V^mol LJ 



which of exactly the same form as (15). A similar calculation in two dimensions reproduces the 
logarithmic dependence in (14). 

A more quantitative comparison to the theories for the VACF tail can be made by examining the 
predictions of the mode-mode coupling theory for the long-time tail, reviewed in detail in Section 
3.2 of Ref. [17]. The relevant formula for the VACF is their Eq. (3.39), which, after integrating 
over the Boltzmann velocity distribution, becomes 



C(t) kLtT 



12vr 3 /9 



/ exp [- (x + v) k 2 t] dk. 

Jk 
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In |47| . the integral over k is performed assuming an infinite system and the time dependence kept 
in order to see the behavior of the tail at long times. If we integrate over t instead, we get 



oc 



Axtaii = / C(t)dt = 3 keT [ k- 2 dk, (25) 
Jt=o 12vr^p(x + v) J k 



which is seen to identical to the integral in Eq. (12) under the assumption that all three directions 

x, y and z are identical (as done in all VACF calculations), 

k B T f kl + kl „ k B T 2 f , _ 2 „ 

{2ttYp (x + v) J k A; 4 (2tt) j p [ X + v) 3 J k 

Note that for a finite system one ought to replace the integrals over k = k ■ 2tt/L with sums over 

k G Z d that exclude n = 0, 

A Xtail = r 2fc f T u3 E fc ' 2 - ( 26 ) 
3p( X + ^3^ o 



In the Molecular Dynamics (MD) literature, the dependence on L _1 in Eq. (24) is considered 
a finite-size effect that ought to be removed in order to extract the bulk (L — > oo) limit of the 
diffusion coefficient [49-51]. A hydrodynamic theory for the finite-size correction, based on the 
Oseen tensor for a finite periodic system, has been developed several times [52] and is confirmed 
numerically in Refs. |51j . This theory focuses on viscous effects only, and we will thus replace v 
with v + D in Eqs. (10,11) in Ref. [51], to obtain 

',2 



6vr (x + v) 



/> max ) + E ^ L ~ 3k ~ 2 ex P (-ji-) 



where f(k) is soniG function. Assuming that A^ max is large, the system-size dependence is captured 



in the last term, which is exactly the same as Eq. (|26|), 

k B T 
6vr (x + v 



1 T 1 

Axmd ~ ( B \ 2 47rL_3fc_2 = Axtaii = A X . 



There are few molecular dynamics studies of the system size dependence of the diffusion coeffi- 
cient in sufficiently large two-dimensional systems. Isobe has performed one of the most extensive 
hard-disk molecular dynamics studies of the hydrodynamic tail of the VACF [53]. We have per- 
formed numerical integration of C{t) using the data of Isobe for hard disks at packing fraction 
4> = 0.18, for square systems sizes containing from = 16 2 to N = 128 2 hard disks. For these 
parameters the statistical accuracy of the data appears sufficient to resolve the asymptotic plateau 
Xeff = lim^oo x(t) °f the time-dependent diffusion coefficient 



X (t) = f C(r)dr, 

Jt=0 
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Figure 8: System-size dependence of the diffusion coefficient for hard disks at volume fraction <j> = 0.f8 (data 
courtesy of Masaharu Isobe [53]). Time and lengths are measured in natural units, or, equivalently, fc^T = I , 
m — 1 and disk diameter a = 1. (Left) The time-dependent diffusion coefficient x(i), as determined from 
a numerical integral of the velocity autocorrelation function, for several system sizes (see legend). (Right) 
The observed increase of the effective diffusion coefficient with system size (symbols), as obtained from the 



large-time limit of x(i). Theoretical predictions are shown for comparison, for both the simple theory (27) 



(dashed line), as well as for the self-consistent theory (29) with the system size N = I6 2 used to define a 



reference length L$. A logarithmic fit with an adjustible slope is also shown (dotted line). 

as illustrated in Fig. [8j Since the aspect ratio of all of Isobe's simulations is fixed at L x /L y = 1, 



Eq. (|20|) suggests that the difference between XeS for a system size of N disks and the smallest 

k B T , N 



system of Nq = 16 disks is 



Xes(N) - Xes{N ) 



hi 



(27) 



8np( X + v) N 

In Fig. [8] we compare this prediction to the numerical integral of Isobe's data, using Enskog kinetic 
theory [53] values for the "bare" transport coefficients x an d v of the hard disk fluid at this density. 
Good agreement is observed, demonstrating that the effect we observe is not an artifact of DSMC 
but rather a generic property of fluids. 



B. Self-Consistent Theory 



The theoretical predictions with which we compared the particle results were based on a leading- 
order perturbative theory |16j that relies on the solution of the linearized equations of fluctuating 
hydrodynamics. A fully nonlinear theory, however, remains illusive. In the context of infinite 
(bulk) systems, a systematic perturbative theory that accounts for corrections of order higher 
than quadratic in the fluctuations has been discussed in Refs. [TT1 [55] . In three dimensions, the 
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conclusion of such studies has been that the higher-order terms do not affect the form of ( 23 ) . In two 
dimensions, several calculations j55H57| and numerical simulations |53} [58] suggest that including 



higher-order terms changes the logarithmic growth in (14). Specifically, it has been predicted that 

-l 

, which 



the self-consistent power-law decay for the VACF is faster than t 1 , C(t) ~ (t\/lnt 
changes the asymptotic growth of Xeff from In L x to \/ln L x . 



In order to obtain a self-consistent form of ( 14 ) , we reconsider the derivation in Section II A 1 



The cause of the diffusion growth from system width L x to L x +dL x is the added contribution to the 



integral in Eq. (13) from wavenumbers in the band \k x \ 6 2ttL~ 1 • [1 — dL x /L x , 1]. If we postulate 
that the concentration fluctuations at wavenumber k x = 2tt/L x evolve with the renormalized 
diffusion coefficient Xeff instead of the bare one, we obtain an ordinary differential equation 
for XcS {Lx), 



dXeS k B T 

dL x (2-K) 2 p(xeS + v)L z 



k x 

l1 



k , 



ky {k x + ky) 



2 dky 



k B T 
4vrp(xeff + v)L. 



r- 1 

^X > 



(28) 



where we have assumed that viscosity does not change substantially with system size (consistent 



with existing molecular dynamics data). Solving the differential equation (28) with the condition 
Xeff (Lq) = x leads to a diffusion enhancement that grows slightly slower than logarithmically, 



Xcs(L x ) « (x + v) 



1 + 



k B T 



In 



L x 



2np{x + v) 2 L z L 



1/2 



v. 



(29) 



Although the self-consistent form (29) still diverges with system length, it is important to observe 



that for any finite system Xeff is well-defined. The predictions of (29) for a hard-disk system at 
packing fraction <p = 0.18 are shown in Fig. [8| We have used the Enskog kinetic theory for the 
viscosity, which is in good agreement with published molecular dynamics data at this density, and 
set pL z to be the mass density in the plane. The self-consistent theory is seen to be in better 



agreement with the molecular dynamics data than the simple theory (14); however, the difference 
between the two is small for the system sizes at which we presently have reliable data for the 
effective diffusion coefficient. 

Repeating the same self-consistent calculation in three dimensions gives the self-consistent form, 

a k B T 



Xeff (L) f« (x + v) 



1 + 







Go 





1/2 



(30) 



which happens to be the solution of the consistency condition 



Xeff = X + 



a k B T 



1 



p[v+ (x + Xcff)/2] VA) L, 
reminiscent of the form obtained in Ref. [T7] except that the arithmetic average of x an d Xeff 
appears in the denominator instead of just Xeff- In three dimensions, the difference between the 
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self-consistent (30) and the simple (15) theories is very small and thus rather difficult to observe 



computationally. 

In both two and three dimensions, the self-consistent predictions ( 29|30 ) will only deviate from 
the simple theory ( 14|15 ) substantially when the diffusion enhancement becomes comparable to 
the bare coefficient, that is, when hydrodynamic effects become comparable to molecular ones. 
The estimates presented in Appendix [A] show that reaching the regime where > X is difficult 
to achieve with particle simulations. In nonlinear fluctuating hydrodynamics finite-volume solvers 
|14j . one has the freedom to choose the various parameters so as to make the effect of advection 
by velocity fluctuations much more prominent, similarly to what is done in Ref. [58J using a two- 
dimensional lattice gas method. Experimentally, two dimensional systems can be realized by using 
thin films, for example, liquid or liquid crystal films. In liquid films, however, velocity fluctuations 
below a certain cutoff wavenumber are suppressed because of the drag by the surrounding fluid, 
and therefore the diffusion enhancement saturates for systems much larger than the cutoff length 
scale [15 1 . 



V. CONCLUSIONS AND FUTURE DIRECTIONS 



The results of our particle simulations confirm that fluctuating hydrodynamics is a powerful 
tool for understanding transport at small scales. Our results conclusively demonstrate that the 
advection by thermal velocity fluctuations affects the mean transport in nonequilibrium finite 
systems and thus the advective nonlinearities, such as the term (5c) (8v) in ought to be retained 
in the equations of fluctuating hydrodynamics. We demonstrated explicitly that the correction to 
the bare or molecular transport coefficients due to the VACF tail [3], hydrodynamic interactions 
with periodic images of a given particle [51] . and the contribution due to thermal fluctuations 
|16l [T8] studied here, are all the same physical phenomenon simply calculated through different 
theoretical approaches, all of which are equivalent because of linearity. The advantage of fluctuating 
hydrodynamics is that it is simple, and it can take into account the proper boundary conditions 
and exact geometry, especially if a numerical SPDE solver is used. Furthermore, other effects such 
as gravity [18], temperature variations [12], or time dependence [3EL5], can easily be included. It 
remains as a future challenge to verify the predictions of fluctuating hydrodynamics for the effect 
of fluctuations on diffusive transport in spatially non-uniform systems, either through particle 
simulations or laboratory experiments [15] . 

Renormalization has often been invoked as a way to fold the contribution from fluctuations 
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into the effective transport coefficients, however, this only works in three dimensions for very large 
systems. In two dimensions, a macroscopic limit does not exist, and in three dimensions there are 
strong finite-size corrections even for systems with dimensions much larger than molecular scales. 
Theoretical modeling of finite systems at the nano or microscale thus requires including nonlin- 
ear hydrodynamic fluctuations. However, a complete nonlinear theory has yet to be developed, 
and requires detailed understanding of the role of large wavenumber cutoffs (regularizations) that 
are necessary to make the SPDEs well-behaved. Furthermore, the proper physical and mathe- 
matical interpretation of other types of nonlinearities in Q and (6|7), notably the dependence of 
the transport coefficients and the stochastic forcing amplitude on the fluctuations, remain to be 
clarified. Future work should also study momentum and heat transfer in steady states, as well as 
time-dependent transport in systems that are far from equilibrium. 



Acknowledgments 

We are grateful to Masaharu Isobe for sharing his hard-disk MD data and helping us analyze it. 
We thank Berni Alder, Doriano Brogioli, Jonathan Goodman and Eric Vanden-Eijnden for informa- 
tive discussions and helpful suggestions on improving this work. This work was supported in part 
by the DOE Applied Mathematics Program of the DOE Office of Advanced Scientific Computing 
Research under the U.S. Department of Energy under contract No. DE-AC02-05CH11231. 



Appendix A: Estimates of the Diffusion Enhancement 

It is instructive to do some scaling analysis of the order of magnitude of Ax in realistic fluid 



systems. Following ( 15 ), the hydrodynamic contribution to the diffusion coefficient for a large three 



dimensional system is estimated as 

Ax ~ -. ^ , (Al) 

PKX + ^jAnol 

For gases, Ax can be estimated by using Chapman-Enskog values for the transport coefficients for a 
hard-sphere gas with molecular collision diameter a « L mo \, specifically, x ~ v ~ (p°' 2 ) 1 \prnki}T ■ 
For liquids, the Schmidt number is large, S c = u/x > 10 2 , and Stokes-Einstein's relation suggests 
that x ~ kfiT/ {pva). For both gases and liquids we get that Ax/x ~ (no -3 ) ~ <p, where n = p/m 
is the number density and 4> is the packing fraction of the particles. We see from this estimate that 
the enhancement due to fluctuations is stronger for dense gases and is strongest for liquids. 
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However, the logarithmic divergence in ( |14[ ) means that the contribution due to hydrodynamic 
fluctuations dominates for sufficiently large (quasi) two-dimensional systems, regardless of the 
density, 

A X k B T L x , 3 . a L x 

~ — — — In ~ (na ) — In — . (A2) 

X PX{X + v)L z L mo i L z a 

A glance at Fig. [6] shows that the enhancement we measured is only a few percent of the kinetic 
theory value (for our DSMC simulations, na 3 = 0.06), and reaching the system width L x where 
A% ~ x is impractical with DSMC simulations at the present. In fact, for the parameters used in 
typical DSMC applications the enhancement of the transport coefficients relative to the Chapman- 
Enskog values is very small. Specifically, assuming an ideal hard-sphere gas collision model and 
taking L mo i = Ax c , for a quasi two-dimensional system of thickness L z = A we obtain the estimate 

X 33ir z (nX 6 ) 

where N\ = n\ 3 is the number of particles per cubic mean free path, and N x = L x / Ax c is the 
number of collision cells along the direction perpendicular to the gradient. In a typical DSMC 
simulation N\ ~ 100, giving A-xlx ~ 0-5 ■ 10~ 3 lnN x , which is less than 0.5% even for N x = 100. 
While using molecular dynamics instead of DSMC allows one to reach larger densities and thus 
enlarge Ax/Xt the regime in which Ax > X is difficult to access even using hard-disk MD [53] (see 
Fig. §. 

In this work we focused on the correlations between velocity and concentration fluctuations. 
Concentration fluctuations also have long ranged self-correlations ~ k~ 4 in the presence of a con- 
centration gradient, see Eq. Q. Even though the enhancement of the concentration fluctuations is 
proportional to the square of the concentration gradient, a two-dimensional calculation [15] similar 
to one presented here [see Eqs. ( 11|13 )] leads to the remarkable result, 

^(-^(-o)g = 12te ,X^. (^)' (Ae) '- 

where Ac = (Vc)L y is the macroscopic concentration variation. Assuming L x ~ L y , we thus obtain 



[c.f. Eq. (A2)] 



ncq K'B-i- 1 3\ °" 



(Ac) 2 PX(X + v)Lz V ' L z 

where na 3 ~ 1 for liquids. We thus see that for systems a few molecules thick, the non-equilibrium 
concentration fluctuations can become comparable to the deterministic variation. In this case we 
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expect that a perturbative approach based on the linearized theory will not be applicable and the 
use of a nonlinear finite-volume solver will be indispensable. 
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